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Abstract. Static earth-tide theory was modihed to include interior ldads expressible as 
spherical harmonics, and elastic moduli were assumed to be functions of radius only. Varia- 
tions of density from this model and the corresponding stress distributions in the crust and 
mantle were calculated to correspond to observed variations in the gravitational field plus the 
surface topography up to 4th-degree spherical harmonics. These solutions were made determi- 
nate by imposing the condition of minimization of strain energy. For the elastic parameters 
derived from seismology, the maximum stress difference obtained from the discrepancy be- 
tween the observed and equilibrium flattenings of the earth was 163 bars. The maximum stress 
difference found for the sum of all other terms was 97 bars in the lower mantle and 300 bars 
in the crust. Displacements were always less than 70 meters. Modifications of the solution 
which take into account finite strain, creep, and viscous deformation are discussed. A model 
consisting of a fluid layer 35 to 400 km deep and a rigidity V 2 seismic in the rest of the mantle 
results in a reduction of the maximum stress difference in the mantle to 54 bars jjg d an in- 
crease of the maximum displacement to 1500 meters. ^ ‘/jS 



Introduction. Knowledge of the earth’s grav- 
ity field has been significantly improved in 
recent years, principally because of the perturba- 
tions of close satellite orbits, but also because 
of the extension of gravimetric coverage and 
the application of large computers to its analy- 
sis. In particular, the long-term perturbations of 
orbits have yielded very accurate determina- 
tions of the low-degree zonal harmonics; the 
most recent summary and analysis were made 
by Kozai [1963]. These zonal harmonics are 
incorporated in the analysis of satellite orbits 
for tesseral harmonics by Kaula [1963]. The 
order of magnitude of satellite results has been 
confirmed by some recent analyses of terrestrial 
gravimetry : least-squares determination of 

harmonic coefficients by Uotila [1962] and au- 
tocovariance analysis determination of degree 
variances by Kaula [1959]. The principal fea- 
tures of variations in the earth’s gravitational 
potential indicated by these studies are: 

1. Zonal harmonics (normalized) up to de- 
gree 4 of the order 10 -6 times the central term. 

2. Zonal harmonics of degrees 5 through 9 
of the order 10' 7 . 

3. Tesseral harmonics up to degree 4 of com- 
parable magnitude to the zonal harmonics. 

4. Negligible correlation with the correspond- 
ing harmonics in the topography for all degrees. 


Stresses in the mantle . Although there re- 
main some statistically questionable aspects of 
all these results except the low-degree zonal 
harmonics, it seems appropriate now to explore 
their implications with regard to stresses in the 
mantle. The subject has been discussed generally 
by O’Keeje [1959], Munk and MacDonald 
[19605], and MacDonald [1963]. Licht [1960] 
has applied the linear viscous model of Vening- 
Meinesz [Heiskanen and V ening-Meinesz , 1958] 
to the explanation of the third-degree zonal 
harmonic. We start at the opposite extreme of 
an elastic model. 

The elastic model is assumed to have elastic 
moduli which are functions of the radius only. 
We seek solutions for variations in density 
which are functions of latitude and longitude 
as well, 8p(r,0,<£), to account for the observed 
variations in the external gravity field and the 
surface load constituted by the topography, 
while at the same time entailing a minimum of 
shear stress in the elastic model. 

The applicable equations in the earth are the 
equations of equilibrium for a continuous me- 
dium: 


0 - 


dW , 



and Poisson’s equation : 


(1) 
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(2) 

where p is density, W is the gravitational poten- 
tial, x i, a*, and Zs are cartesian coordinates, 
is the stress tensor, and G is the gravitational 
constant. 

The potential can be considered as consisting 
of a potential existing in a reference state plus 
a perturbation : 

W = Wo + ¥ (3) 

The corresponding equation for the density is 

p — Po ~ J diii — p 0 A + bp (3a) 

where A is the dilation and Bp is a density per- 
turbation from outside the problem: a chemical, 
structural, or thermal inhomogeneity. In (3a), 
the rule of summation over the repeated sub- 
script applies, and it will apply in all subsequent 
equations. 

We make the assumptions that: (1) the ref- 
erence state is one of fluid equilibrium, (2) the 
stress-strain relationship is elastic, and (3) the 
displacements u< = x t — x oi are small enough 
so that only linear terms must be considered. 

These assumptions enable us to write the con- 
stitutive equation as 

Pii = SijXA + M (& + g) (4) 

X = k - 2/i/3, (4a) 

where 8*, is the Kronecker delta, /x is the modu- 
lus of rigidity, A is the Lame constant, and k 
is the bulk modulus. 

The problem thus far is identical with that of 
earth tides, which has been treated by Takeuchi 
[1950], Molodenskii [1953], and Alterman et al. 
[1959], with one exception: the term Bp in the 
expression for the density, (3a). Substituting 
the expressions from (3) and (3a) into (1) and 
(2) and neglecting the products of small quan- 
tities, we obtain terms additional to those of the 
earth tide problem of 

bp dWo/dXi (la) 

on the right of (1) and 

“47 r G bp (2a) 

on the right of (2). 

If the hydrostatic reference state is subtracted 
out, and the equations converted to spherical 


coordinates, as in equations 7 through 10 of 
Alterman et al. [1959], the density perturbation 
will appear only in the radial equation of equi- 
librium (7) as 

“ &P Qo (1 & ) 

where g 0 is the negative of the radial gradient 
of the reference potential. Following Alterman 
et al. [1959] further, we express the displace- 
ment as a spheroidal vector spherical harmonic : 


= U(r) 


0 
0 

(A(0, <p)i 


+ V(r) 


dS n (e, <p)/de 

[dS n (e f <p)/d<p]/sm 0 

0 


(5) 


where S n ( 6 , <p) is a surface spherical harmonic, 
and we express the perturbation of the potential 
as 


* = P(r)S n (S t <p) (6) 

Furthermore, we express Bp as 

5 P - D(r)S n (6,<p) (7) 

Substituting (5) and (6) into equations 7 
through 10 of Alterman et al. [1959], and sub- 
stituting (7) for our extra terms (16) and (2a) 
we get the additional term — Dg 0 on the left 
side of their equation 23 and the additional term 
— 4 t tGD on the right of their equation 25. 

Finally, following their conversion of variables, 
we get the six first-order equations: 


Vi 




2 \yi 

(X + 2 p)r 


+ 


V2 


+ 


Xn(n + 1) 


(X + 2 p) (X + 2fi)r yz 

' [ s | 4m(3X + *)> 

V 2 = [- 4 po! 7 or+ (x+ Y m) -J ? 


4/i 


n(n 


y 2 + n(n + l)p c g 0 r 


(X + 2fi)r 1 

_ 2^(3X + 2ju)n(n + l) ~]y 3 
(X + 2 P ) Jr’ 

, n(n + 1) , ^ 

H ~ Vt — PoVt + Dg 0 
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Ve 


V* 


_yx + vi + y± 

r r p. 

[~ 2m(3X + 2p)~\yi 

I*" ~ (X+2„) J? 

X J/2 , 2 fi 

(X + 2{j) r (X + 2/x) 

• [X(2n 2 + 2n - 1) 

+ 2 n(n 2 + n — 1)] 

r r 


PoV 5 

r 


( 8 ) 


?// = 4:irGp 0 yi + ye 
, 4irGp 0 n(n + 1 ) 

Ve = 2/3 

r 

, n(n +1) 2y 6 A ^ 

H 2 $/s 47 rCrZ) 

r r 

where the y <’ s are the radial factors, respectively, 
of the radial displacement, the radial stress, the 
tangential displacement, the tangential stress, 
the potential perturbation, and the gradient of 
the potential perturbation less the radial dis- 
placement contribution thereto. Primes denote 
derivatives with respect to r. 

We abbreviate (8) as 


Vi' = QaVi + WiD (9) 

At least one Q i} is of the order of 1/r for all j, 
so all y / s must be zero at the origin. Further- 
more, some Q ilf Q i9) and Q iS are 0(l/r*), so 
Vi'(O), y 3 '(0), and y/{ 0) must all be zero, leav- 
ing three constants of integration, y 2 '( 0), y/( 0), 
and y a '( 0), for a solid core. For a fluid core, 
there are only two independent equations (cf. 
Longman [1963]): 




y e 




n(n + 1) 
2 

r 




( 10 ) 


For a fluid core with radius c, we took as con- 
stants of integration i/«'( 0), t/i(c), and t/ s (c). 
The conditions at the surface r = a with a 
layer of surface density <r n are 


y 2 (a) = — <7o(a)cr n | 
2/4 (a) = 0 

[(n + 1)/ a] 2/5(0) + 2/s(a) = 47rG(r„ 


If we consider the external gravitational field 
as known, with coefficient C„ for the spherical 
harmonic component of the potential, there is 
a fourth surface condition, 

2 h(d) - C n (11') 

If it is also assumed that D(r) is known 
throughout, the problem is overdetermined. The 
simplest assumption is that D(r) is known ex- 
cept for a constant multiplier k, which becomes 
a fourth unknown in addition to the three con- 
stants of integration. If <r„ is zero and C n is the 
exterior potential due to D(r) in a perfectly 
rigid earth, the load Love number k n " for the 
internal load D(r) comparable to the cal- 
culated for surface loads by Takeuchi et al. 
[1962], and Longman [1963] is 

k n " = (1 - k)/k (12) 

k, like k", is thus a measure of the response of 
the earth to the internal load Bp] the exterior 
potential perturbation C„ will be a combination 
of the potential due to Bp itself and that due to 
the deformation of the earth by the load Bp. 

If D(r) is assumed to involve two or more 
unknown parameters, further assumptions are 
required in order to make the problem de- 
terminate. The most logical assumption is that 
the shear strain energy is minimized; i.e., 

J ju(r)ei/e { / d V = minimum (13) 

where ju is the rigidity, the integration is over 
the volume of the earth, and e</ is the deviatoric 
strain tensor [ Jeffreys , 1959 p. 12] 


&i j j 3^1 j &kk (id) 

To connect e if with the y/ s, we use the equa- 
tions for the strain in spherical coordinates 
[Love, 1927, p. 56], (applying a factor to the 
off-diagonal components of strain to be con- 
sistent with tensor convention) and eliminate 
dyjdr and dy*/dr by (8) 


V\ + ~ 

r r 


2 d 2 S n . 


ar 


y% 


e v<p 



2 

r sin 6 




+ cos 6 


dd ) 


( 11 ) 


2/3 
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e rr = - 


_4X_ Sn, 
X + 2 ii r 


V\ + 


2n(n + 1)X 
(X + 2y) 


2 S n . 


$nm I 

•T" 2/3 + X + 2m 


2/2 




€ vr 


e r e 


2 ( d*s nm _ 


r sm 0 \d0d<p 

1 


cot 0 


dS n „ 

dip 


2/2 


fi sin 0 d<p 

i as,. 

At dO 


2/4 


2/4 


We express (15) in matrix form 

e(0, <P, r) = M(0, ¥>)y(r) 


(15) 


(16) 


The shear strain energy per unit volume can 
then be expressed as 


«(*,*, r) = mW y r M T JMy (17) 

where the r superscript denotes the transpose 
and 


! -a -i o o o' 

-i § -i o o o 

J= -1 -i ! o o o 

0 0 0 2 0 0 

0 0 0 0 2 0 

. 0 0 0 0 0 2, 

The Runge-Kutta or other numerical method 
of solution of (9) (see, e.g Royal [1955]) can 
be adapted to express the variables ?/ { at level 
k in terms of the y t at level k — 1 and the 
density function D(r) at levels k — 1, k — 4, 
and k: 


y k = T k y k - X + w k D k + w 2 *D*_ 1/2 + w 3 *D*-i 

(18) 

If the density function D(r) is expressed as a 
function of a few parameters d, (18) can be 
simplified to the form 


y* = T^y*., + W t d (19) 

By repeated substitution of the k — 1, k — 2, 
etc., equations into (19) we obtain 


( 20 ) 


where c is a vector of the constants of integration, 
V>( 0), Vi(c), and y 3 (c), 

Uk = fi T, 

J-l 

and 

F* = Z n Tj-W,-, + W* 

t-2 j-t 

and Ti is appropriately modified to fit the in- 
itial conditions. 

The surface conditions (equations 11 and IT) 
for m levels of integration are expressible as 

a = Cy m = CU m c + CV w d (21) 

The system of equations (17), (20), and (21), 
along with the condition expressed by (13), is 
mathematically identical to that of generalized 
least squares (see, e.g., Arley and Buck [1950, 
pp. 196-198]) with the vector of corrections to 
observations 


x = (c^dy (22) 

condition equation coefficient matrix 

F - (CU m • CV m ) (23) 

and covariance matrix 

w = EmW(u t -v t y 

k 

•( / M T JM d»)(U 4 - • • V*) (24) 

where the integration is over the surface of the 
sphere. The solution of the system of equations 
(21) through (24), subject to (13), by the 
method of Lagrangian multipliers is 

x = WF T (FWF T ) -1 a (25) 

Solutions were made for density distributions 
corresponding to spherical harmonics of the ex- 
ternal gravitational field up to the harmonic 
7 M given by Kaula [1963] and the Gutenberg 
model density and elasticity parameters as given 
by Takeuchi et al. [1959]. The surface layer 
coefficients were derived from the harmonic 
analysis of the topography by G. J. Bruins, as 
described by Vening-Meinesz [1959]. In these 
solutions the density parameters for the mantle 
and the crust were kept separate. Two different 
types of parameters d were used to represent 


y* — c + V^d 
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TABLE 1. Mantle Model Corresponding to the External Gravitational Field 
Minimum shear strain energy, Gutenberg model of elastic parameters. 


Coefficient 

External 
Potential, 
earth units 

Topography 
Surface 
Density, 
equivalent 
earth units 

Crustal 

Density 

Anomaly, 

g/em* 

Maximum 

Mantle 

Density 

Anomaly, 

g/em* 

Maximum Mantle 
Stress Difference 
at Various Radii 

Total 

Shear 

Strain 

Energy,* 

ergs 

Ap, R, 

bars earth radii 

AC 20 

5.0 X 10-« 

4.65 X 10" 6 

0.012 

5.7 X 10- 

4 163 

0.56 

22.4 X10 28 

C22 

1.84 

-3.9 

0.021 

-1.2 

70 

0.56 

4.6 

>§22 

-1.71 

-0.34 

-0.002 

-1.3 

65 

0.56 

2.4 

CsQ 

0.98 

-2.50 

0.014 

-1.4 

32 

0.56 

2.4 

Cii 

1.77 

-1.52 

0.012 

-2.2 

60 

0.56 

5.4 

S n 

-0.11 

1.16 

-0.006 

1.2 

5 

0.78 

0.1 

Cj2 

0.34 

-4.45 

0.022 

-2.9 

15 

0.56 

2.3 

jSj? 

0.08 

3.94 

-0.018 

2.7 

3 

0.56 

1.6 

C33 

-0.31 

0.72 

-0.004 

0.3 

11 

0.56 

0.2 

Ssz 

0.74 

4.47 

-0.019 

3.3 

55 

0.56 

2.8 

C40 

-0.41 

2.68 

-0.014 

2.1 

18 

0.78 

0.6 

C 41 

-0.21 

-1.68 

0.007 

-1.2 

8 

0.87 

0.4 

&41 

0.46 

-2.46 

0.013 

-1.8 

19 

0.77 

1.2 

Cj2 

-0.03 

-4.0 

0.018 

-2.9 

9 

0.98 

1.6 

Sn 

0.32 

0.6 

-0.002 

0.5 

12 

0.56 

0.1 

Cn 

0.50 

3.0 

-0.012 

2.2 

15 

0.86 

1.4 

*§43 

0.16 

-1.8 

0.009 

-1.3 

7 

0.98 

0.4 

C44 

-0.24 

-0.19 

0.000 

-0.2 

9 

0.56 

0,2 

Su 

0.55 

4.25 

-0.018 

3.1 

20 

0.56 

2.5 


* Neglecting interactions between different harmonics. 


the density anomalies in the mantle: layers and 
polynomial coefficients. About the same answers 
were obtained by the two methods, but those 
from the polynomial coefficients varied more 
smoothly. In general, no appreciable reduction 
in the summation of strain energy was obtained 
by using more than four parameters to repre- 
sent the density variations for each harmonic 
in the mantle. One parameter was used to 
represent the density anomaly in the crust. 

As a test of the program, tidal and surface- 
load Love numbers were also calculated by fix- 
ing D(r) = 0 and appropriately modifying the 
surface conditions. The answers agreed more 
closely with those of Longman [1963] than 
with those of Takeuchi et al. [1962]. 

The results for the polynomial density-varia- 
tion solution are summarized in Table 1. The 
potential and surface density coefficients are in 
units such that the radius of the earth, the 
mass of the earth, and the gravitational con- 
stant are all unity, and all coefficients apply to 
normalized spherical harmonics S nm so that the 
integral of S nm * over the unit sphere is 47r. The 
AC* is the discrepancy between the observed 


value and that corresponding to fluid equilibrium 
[O'Keefe, 1959] . The maximum shear stress ob- 
tained for AC* is consistent with the 100 bars 
strength estimated by Munk and MacDonald 
[1960a, p. 280]. 

Figures 1 and 2 sum up the effects of all coef- 
ficients except AC 20 in the form of maps of the 
maximum stress difference and the radial com- 
ponent of displacement at three selected levels 
within the mantle. The displacements have, in 
general, a negative correlation with the varia- 
tions in the external field [Kaula, 1963, Figure 
1]. The stress differences, which are quadratic 
functions of the displacement gradients, do not 
show such a clear pattern, although maximums 
generally occur in areas where the correlation 
between the topography and the gravitational 
field is most negative. The appreciable varia- 
bility of the stress differences suggests that a 
significantly different solution might be obtained 
by applying a yield stress limit. Also calculated 
were strain energy densities ; the maximum shear 
strain energy density found in the mantle was 
1100 ergs/cm*. 

As might be expected under the criterion of 
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Fig. 1 . Principal stress differences in bars, rigid- 
ity and bulk modulus assumed from seismology 
(Gutenberg model). 
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0.6 

EARTH 

RADII 


R= 

0.83 
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R= 

0.98 
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Fig. 2. Radial displacements in meters, rigidity 
and bulk modulus assumed from seismology (Gu- 
tenberg model). 
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Fig. 3. Principal stress differences in bars, bulk 
modulus assumed from seismology throughout ; 
fluid layer 35 to 400 km deep ; and Y 2 seismic 
rigidity 400 to 2900 km deep. 
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Fig. 4. Radial displacements in meters, bulk 
modulus assumed from seismology throughout; 
fluid layer 35 to 400 km deep; and V 2 seismic 
rigidity 400 to 2900 km deep. 
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strain energy minimization, for most terms the 
density anomaly in the crust accomplished al- 
most complete isostatic compensation of the sur- 
face topography. However, in satisfying the 
condition of the known external gravitational 
field the compensation is never exact : generally, 
for terms in which the surface layer has a sign 
opposite to that of the gravitational field term, 
there is an overcompensation, and if it has the 
same sign, there is an undercompensation. The 
maximum density anomaly in the mantle nearly 
always fell immediately below the crust. 

A manifestly oversimplified physical model has 
two possible values: to give a notion of the 
quantities involved and to suggest how a more 
realistic (and hence probably more complicated) 
model might be developed. The stresses and 
energies found for the elastic model are prob- 
ably a crude approximation to those in the 
actual earth, but the deformations obviously 
fall far short of explaining what occurs. We 
therefore want to examine in turn the possible 
modifications of the model in the categories of 

(1) finite strain elasticity: the strain is con- 
stant, but the stress is significantly affected by 
higher-order terms than those in equation 4; 

(2) creep: the strain rate is appreciable, but 
the stress remains a function of only the strain; 
and (3) viscous deformation: the stress is a 
function of the strain rate. 

Finite strain elasticity. Finite strain is sig- 
nificant within the earth in two respects: first, 
the density and elastic moduli for small super- 
imposed loads, such as seismic waves, must vary 
with pressure and temperature in a way con- 
sistent with actual materials; second, the higher- 
order terms may contribute appreciably to the 
stress: i.e., products of displacement gradients 
may be significant in calculating strain, and 
products of strain components may be significant 
in calculating stress. 

The first aspect of the finite strain theory 
has been applied most extensively by Birch 
[1952, 1961] to the problem of deducing com- 
position or phase changes in the mantle from 
seismic velocities. It is of concern in our pres- 
ent problem only to the extent that the density 
and elasticity parameters assumed should be 
materially plausible. 

The second aspect of finite strain theory seems 
of dubious applicability to the earth's mantle, 
since the strains involved in the minimum energy 


solution described above are less than 10"\ If 
higher-order terms were taken into account, the 
initial state of stress would no longer be simply 
additive; hence the solution according to strain 
energy minimization or other criterion would 
be a function of the initial state of stress. This 
problem was investigated by Jeffreys [1943], 
who assumed an implicit nonlinear constitutive 
equation in the form of a stress matrix with 
one higher-order term of the form xa The 
model was a homogeneous shell with density 
layers at the upper and lower boundaries. The 
radial factors of the displacements were ex- 
pressed in terms of four functions of the form 
A, -f B,f. The eight parameters A s and B, 
were adjusted to minimize the strain energy or 
stress difference while satisfying the boundary 
conditions separately for each harmonic. The 
reduction in stress obtained was moderate, # even 
though interactions between different harmonics 
were neglected and it was not made clear 
whether the constitutive equation implied by 
the solution was realistic. 

In our formulation of the problem, it seems 
very probable that considerations of nonelas- 
ticity are of much greater concern that nonlinear 
elasticity terms. 

Creep. There have been several recent dis- 
cussions of the geophysical application of ex- 
perimental evidence and theoretical models of 
creep: Orowan [1960], Griggs and Handin 
[1960], Jeffreys [1958, 1959], Jeffreys and 
Crampin [1960], MacDonald [1961], Lomnitz 
[1962], Weertman [1962], Scheidegger [1963], 
Stacey [1963], and Donath and Fail l [1963]. 
There appear to be appreciable differences of 
terminology, interpretation, and opinion in these 
discussions, so we attempt to summarize the 
principal conclusions in order to decide how our 
mathematical model should be modified. 

1. Transient creep, or elastic afterworking. 
Consideration of the phase lags and dampings 
in the response of the earth to periodic dis- 
turbances leads to creep models of the form 
[ Jeffreys , 1958; Jeffreys and Crampin, 1960; 
MacDonald , 1961 ; Lomnitz , 1962] 

e(<) = PM _1 [1 + 'KO] (26) 
for the strain e at a time t > 0 due to a con- 
stant stress p applied at t = 0. 'Creep functions' 
'Sr(t) which appear to fit phenomena having 
periods up to the 430 days of the free nutation 



4974 


WILLIAM M. KAULA 


are 




\f/(t) = q log (1 + at) 

(27) 

and 




'P(i) = [(1 + at) a — 1 \q/a 

(28) 


The time scale of these phenomena falls far 
short of that in our problem, but Jeffreys and 
Cr ampin [1960] find that their numerical values 
in (28) would permit an appreciable part of a 
second-degree harmonic to survive for more 
than 10 8 years. MacDonald [1961] objects to 
the rules (27) and (28) because they imply an 
infinite population of relaxation times in a linear 
superposition model and because they do not 
lead to frequency independence of dissipation, 
as observed for frequencies above 1 cps, but 
these considerations have no apparent bearing 
on tjie long term problem. 

2. Delayed elasticity. Glasses at low tempera- 
ture can be regarded as elastic, but their full 
response to loads is appreciably delayed. If the 
irregularities of grain boundaries, etc., are such 
that the material of the mantle can be regarded 
as amorphous, there may be a significant 
amount of delayed elasticity. Orowan [1960] 
suggests that delayed elasticity may account 
for the post-glacial rise. In this case, elasticity 
would have to prevail to a considerable depth in 
the mantle, rather than only in the crust, as 
considered by Vening-Meinesz [ Heiskanen and 
Vening-Meinesz, 1958]. 

3. Steady-state creep, or elasto viscosity. 
Crystalline substances at low shear stresses 
above a 'creep strength’ deform very slowly. 
This steady-state creep can usually be ex- 
pressed as [Weertman, 1962; Stacey, 1963] 

e = C(p) exp ( — Q/kT) (29) 

Where e is the strain rate, Q is an activation 
energy, k is the Boltzmann constant, T is the 
absolute temperature, and C is an exponential or, 
near melting, a power of the shear stress p. 
Jeffreys [1958] considers steady-state creep to be 
of no importance in the earth; Stacey [1963] 
suggests that it is the dominant rheological mode 
in the mantle. 

4. Uniform flow. If the shear stress in a 
crystalline substance exceeds a certain yield 
point, rapid plastic deformation or uniform 
flow takes place. At confining pressures above 
1000 bars, the zone of failure widens, but cohe- 


sion is retained after failure [Griggs and Han- 
din, 1960; Donath and Faill, 1963]. 

Instability of creep in metals has been ob- 
served to result in the rapid propagation of 
intense shear bands, or faulting without fracture 
[Orowan, I960]. This 'shear melting’ is sugges- 
tive of a mechanism for earthquakes or magma 
formation. 

Estimates of the creep strength or yield point 
at the temperatures and pressures prevailing in 
the mantle seem to be rather speculative. Orowan 
[1960] states that the creep strengths of crys- 
talline materials close to the melting point are 
of the order of 100 or 10 bars, or even less. 
Stacey [1963], on the other hand, obtains creep 
strengths (defined as the 'stress at which creep 
reaches [a] just observable value’) in excess of 
500 bars, mainly, it seems, from the observa- 
tion of Scheidegger [1963, p. 159] that depar- 
tures from isostasy which result in stresses 
smaller than about 4000 bars can exist 'almost 
indefinitely.’ Weertman [1962] does not give a 
creep strength, but he does say that a steady- 
state creep law of the form of (29) breaks down 
at stresses above 100 to 1000 bars, which pre- 
sumably defines the yield point. Griggs and 
Handin [1960] make no sharp distinction be- 
tween creep and uniform flow but suggest that 
shear melting could occur in the mantle at 
shear stresses as low as 100 bars. 

We seem to have authority for a wide range 
of hypotheses as to creep strengths and yield 
points in the mantle. In particular, in compar- 
ing the stress differences in Figure 1 with the 
values suggested for the strengths, it appears 
entirely possible that the mantle can be as- 
sumed to be in a state in which creep is taking 
place, but in which stress is a function of the 
strain, not of the strain rate. Such a state would 
occur if the density irregularities Sp were not 
removed by the creep in the 10 8 or 10 9 years 
time scale of our problem. Taking the situation 
at any instant t = 0 as a reference state, we 
find that the rigidity p, in (4) would be replaced 
by a pseudo-rigidity 

nM = n/[i + m (30) 

where 

i p(t) = [ A(t — r)p(r) dr > 0 
do 

which is formally similar to the transient creep 
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TABLE 2. Mantle Model Corresponding to the External Gravitational Field 
Minimum shear strain energy, fluid layer 35 to 400 km deep, \ Gutenberg model rigidity in lower mantle. 


Coefficient 

External 
Potential, 
earth units 

Topography 
Surface 
Density, 
equivalent 
earth units 

Crustal 

Density 

Anomaly, 

g/em 3 

Maximum 

Mantle 

Density 

Anomaly, 

g/cm 3 

Maximum Mantle 
Stress Difference 
at Various Radii 

Total 

Shear 

Strain 

Energy,* 

ergs 

Ap, 

bars 

R, 

earth radii 

A (?20 

5.0 X 10~ 8 

4.65 X 10-« 

0.0146 

- 2.0 x 10- 4 

31 

0.56 

1.84 X 10 !8 

(?22 

1.84 

-3.9 

0.0019 

-1.0 

13 

0.56 

0.19 

$22 

-1.71 

-0.34 

-0.0042 

0.7 

12 

0.56 

0.19 

<?30 

0.98 

-2.50 

0.0034 

-3.5 

22 

0.937 

0.40 

C« 

1.77 

-1.52 

0.0080 

- 6.1 

35 

0.937 

1.39 

$31 

- 0.11 

1.16 

0.0002 

0.4 

2 

0.937 

0.05 

Cj 2 

0.34 

-4.45 

0.0010 

-1.4 

7 

0.937 

0.06 

4 82 

O.US 

3.94 

0.0028 

0.03 

1 

0.56 

0.03 

Cil 

-0.31 

0.72 

- 0.0011 

1.1 

6 

0.937 

0.04 

Su 

0.74 

4.47 

0.0065 

- 2.3 

13 

0.937 

0.32 


-0.41 

2.68 

-0.0013 

2.7 

16 

0.937 

0.14 

Cji 

- 0.21 

-1.68 

-0.0024 

1.2 

7 

0.937 

0.05 

Si l 

0.46 

-2.46 

0.0017 

-3.0 

16 

0.937 

0.08 

c i2 

-0.03 

-4.0 

-0.0025 

-0.09 

5 

0.82 

0.03 

S 42 

0.32 

0.6 

0.0025 

-1.9 

9 

0.937 

0.10 

Cj 3 

0.50 

3.0 

0.0052 

-2.9 

13 

0.937 

0.26 

*§4| 

0.16 

- 1.8 

0.0006 

- 1.1 

5 

0.937 

0.02 

C 44 

-0.24 

-0.19 

-0.0018 

1.5 

7 

0.937 

0.05 

Su 

0.55 

4.25 

0.0062 

-3.1 

16 

0.937 

0.34 


* Neglecting interactions between different harmonics. 


function of MacDonald [1961]. Equation 30 
makes ii p (t) a function of latitude and longitude, 
which in turn would not permit development of 
the solution in the form of (8). To see how the 
stresses and displacements might be modified 
in an extreme case, however, a solution of the 
system of equations 21 to 24 subject to (13) 
was made for a model in which jx p (t) was zero 
in the low-velocity zone 35 to 400 km deep, 
and equal to the Gutenberg model rigidity 
throughout the rest of the mantle. This solution 
is summarized in Table 2. The principal result 
was that the incompressibility of the fluid layer, 
in addition to that of the core, reduced the 
shear stresses required by. supporting much of the 
load. Specifically, there was (1) an increase in 
the radial displacements of the crust by a factor 
of about 30, (2) a decrease in the shear stresses 
in the crust to about (3) a reduction of the 
stresses in the mantle to about £ for the second 
degree, \ for the third degree, and £ for the 
fourth degree, (4) a concentration of stresses 
in the mantle near the solid-liquid boundaries, 
and (5) a density-anomaly maximum at the 
top of the elastic part of the mantle. The re- 
duction of the ACao stress to a value comparable 


to that for C S1 suggests that reduction of rigidity 
similar to that of the hypothesized model could 
take place in the 10 7 years indicated by the lag 
in adjustment of the rotational bulge. The rela- 
tively small reduction in stresses for the fourth- 
degree terms suggests that higher-degree varia- 
tions could not be reasonably supported by such 
a model, and it is consistent with the indications 
from satellite zonal harmonics and auto covari- 
ance analysis of gravimetry that there is an 
appreciable drop-off in the magnitude of varia- 
tions of the gravitational field above the fourth 
degree. 

Viscous deformation . We may ask two ques- 
tions: Is large-scale convective motion now 
taking place in the mantle? If convection is 
taking place, is the stress which supports the 
density irregularities primarily a function of the 
strain rate — i.e., is it viscous (Newtonian or 
otherwise)? Answers to the first question appear 
to depend on the background of the answerer, 
the two extremes being the geothermists [Lubi- 
mova, 1960; MacDonald , 1963], who consider that 
convection must be limited in time and place 
(of the nature of shear melting) in order that 
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the observed heat flow not be greatly exceeded, 
and the paleomagnetists [ Runcorn , 1962; Stacey, 
1963], who need to explain evidence of conti- 
nental drift of the order of 3 cm/year. Answers 
to the second question are even vaguer because 
of the lack of experimental data, the uncertainty 
as to the appropriate rheological theory, and 
the lack of an adequate mathematical solution 
for convection with significant variation of 
parameters and a finite yield stress, as discussed 
by MacDonald [1963]. If it is assumed that stress 
is a function of strain rate d {i with constitutive 
equation p {i = pd kk bu + 2 tj d ijf where rj is the 
viscosity, a solution for the displacement rate 
field analogous to (4) to (25) for the displacement 
field would be somewhat more complicated, 
since the external gravitational field would still 
be a function of the displacement field. We should 
expect, however, that the order of magnitude of 
the rates v could be estimated from the displace- 
ments u by 

v uy/rf (31) 

The pseudo-viscosity rj most commonly 
quoted is 10 22 poises, deduced by Vening-Meinesz 
[ Heiskanen and Vening-Meinesz , 1958, pp. 365- 
370] and others from a decay time of about 
5300 years and a linear dimension of about 1200 
km for postglacial uplift. The Newtonian 
viscous model requires that the decay time vary 
inversely as the diameter of the load; however, 
Crittenden [1963] recently obtained a decay 
time of about 4000 years for the 180-km-diame- 
ter Bonneville Lake area uplift, or a pseudo- 
viscosity r) of 10 21 poises. 

Taking y = 10 12 , u — 3 X 10 3 , and r? = 10 2S 
cgs in (31), we obtain velocities of the order of 
10 cm/year. Licht [1960] objected to such a high 
rate because it required improbably high effi- 
ciency of heat transport from core to surface; 
however, his discussion seems to be based on the 
implicit assumption that the radiogenic heat in 
the mantle is less than 10 -2 of the 1.6 ergs/g/yr 
estimated for chondritic composition [ MacDonald , 
1959]. Runcorn [1962] also appears to neglect 
radiogenic heating in the mantle in showing that 
convective velocities of the order of 10 em/yr 
are obtained with temperature differences of only 
0.2 degree centigrade, thus yielding a surface 
heat flow much less than observed. A more 
fundamental weakness of mantle-wide convection 
models is the assumption that the pseudo- 


viscosity of 10 22 poises can be applied to organized 
convection with a characteristic length of 2900 
km; Crittenden's result emphasizes the danger 
of such an extrapolation. 

To reconcile the gravitational and thermal 
evidence with a Newtonian viscous model ap- 
pears to require either a mantle-wide system 
with viscosity coefficient in excess of 10 23 poises 
or a system with a characteristic length of 180 
km or less and a viscosity of 10 21 poises or less. 
The first model leaves unexplained many sur- 
face evidences of large-scale motion ; the second 
model has not been worked out. Therefore, 
further development of convective models seems 
worth while, as well as critical re-examination 
of the evidences of continental drift [ Munk and 
MacDonald , 1960a, pp. 251-262, 282-285] 
and study of the problem of distribution of 
radiogenic heating. 

Conclusions . The indications we have found 
as to the state of stress in the earth's interior 
from the low-degree, or long- wave, variations 
in the external gravitational field are minimal in 
two respects. First, the condition of strain en- 
ergy minimization has been imposed ; the actual 
density distribution could conceivably be quite 
different. Second, superimposed on the low-de- 
gree variations may be higher-degree variations 
of sufficient magnitude to increase appreciably 
the stress differences above those here calcu- 
lated. If the scale of significant change in the 
earth is appreciably smaller than the wave- 
lengths of the low-degree harmonics — as sug- 
gested by the narrowness of the belts of surface 
manifestation of tectonic activity and the rela- 
tively shallow origin of seismic activity — then 
the long-wave variations in the gravity field 
supply evidence of only the passive background 
for the processes currently important. 

The development of models more realistic 
than the elastic cannot be purely mechanical be- 
cause of the rheological uncertainties. To pro- 
vide some limitation to the possible solutions, we 
should include energy flows and distribution of 
heat sources in the problem. The desirability of 
incorporating the thermal aspects, as well as 
shorter-wave variations in the gravitational 
field, suggests a more statistical approach to 
keep the problem manageable. The inputs would 
be the spectrums of gravity, topography, heat 
flow, etc., variations in the form of degree 
variances 
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<?n = L (.CnJ + S nm 2 ) 

m— 0 

together with some measure of the correlation 
between the different observed quantities. The 
solutions sought would be the speetrums of 
density, displacement, displacement rate, heat 
sources, etc., at various levels within the earth 
under specified conditions such as the strain 
energy minimization. 
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